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0 Introduction 

We will consider the second-order linear differential equation 

(py'Y + qy = (Ain + A 2 r 2 h -h A d r d )y ( 1 ) 

on a real interval X\ < x < X 2 , where p, q, n,..., r<j are given functions and Ai,..., A d are 
unknown parameters. Let some appropriate boundary conditions be imposed at x\ and X 2 - 
Then a spectral problem consists of determining the subset (Ai,..., Xd) C R d (or C d ) for which 
there exists a solution y of © which satisfies those boundary conditions. While there is a 
vast literature on spectral theory for general differential equations and on numerical methods 
specifically developed for m for d = 1- indeed, the expression “spectral problem” commonly 
implies a single Ai — there is considerably less available concerning several parameters. One 
may find qualitative results on this subject in m ana m Ea sa. For d = 2 some properties of 
the eigencurves are set forth in [31 chapter 6]. 

An approach for solving spectral problems for d = 1 was presented in 1201 1231 which produces 
two explicit power series in the variable Ai with coefficients which are functions on [x\, £ 2 ]. 
These series represent two functions y = u\{x), y = U 2 {x) parametrized by Ai which are linearly 
independent solutions of ©■ There are similar power series for the derivatives uRa;), u' 2 {x). 
By evaluating these series with x at the endpoints x±, X 2 we obtain power series in Ai for the 
boundary values, which upon substitution in the boundary conditions produce a “characteristic 
function” whose zeroes are the eigenvalues of the spectral problem. (It is not necessary for the 
boundary conditions to be linear for this procedure to apply.) These series representations have 
applications beyond spectral problems; for example they provide an effective method for solving 
initial value problems. 
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Since its appearance in 2008, consequences of this SPPS (spectral parameter power series) 
representation have been investigated in many directions. These include completeness properties 
of the “formal powers” used to define the coefficients of the power series mm, relationship to 
transmutation operators, Darboux and other transformations, and Goursat problems m m m 
[28] ; extension to other number systems (quaternions, etc) 0 OS GET] and equations of higher order 
m-, relaxation of regularity conditions on the coefficients of the differential equation 01111- 
Further, there have appeared numerous applications to problems in physics and engineering 
[CL EB iIEj \m 122] as well as in complex analysis [5] [Ml . 

In dealing with physics or engineering models which involve a Sturm-Liouville problem 
containing several eigenvalues A;, it is common practice to fix all but one of them, and then 
solve the spectral problem for the remaining one. This appears to be due to the difficulties 
of existing methods of handling more than one spectral parameter. In this paper we work 
out the SPPS coefficients corresponding to © for arbitrary d > 1. This permits treating the 
spectral parameters in unified way. We give some numerical examples with d — 2, and then an 
application to a problem of transmittance of an electromagnetic wave through an inhomogeneous 
layer, in which the two parameters correspond to physical characteristics of the phenomenon. 


1 Formal powers 


The Sturm-Liouville linear differential expression on the left-hand side of © will be denoted by 


Ly = ( py'Y + qy- 


( 2 ) 


Throughout this paper p, q, rq,..., will denote real or complex valued functions on the closed 
interval [aq, aq]. In this section we are interested in describing the procedure for constructing 
the SPPS representation of solutions, while questions of convergence and regularity will be 
deferred to the next section. A basepoint xq is fixed in [aq, X 2 ] ■ For convenience we will use the 
notation 



to mean 



f(s)ds 


for any function / under consideration, inasmuch as we will have no use for other limits of 
integration. In the following we will set up a notation for describing sums of finitely nested 
integrals of the form 
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1.1 Construction of 


For simplicity of handling the indices, we will begin with the form ([3| which produces the family 
of functions called X^' 1 (the notation follows that generally used in the SPPS literature). Here 
j = (ji.... ,jd) is a multiindex with integral entries. It has d predecessors given by 

J &i = (jl) • ■ ■ ) ji— 1) ji lj ji+11 ■ ■ ■ tjd) 

where 5i = (0,..., 0,1, 0,..., 0) is the i-th standard basis vector. We will say that j is an 
admissible multiindex when at most one of its entries ji is odd. An admissible j is called even 
or odd according to the parity of 

d 

IjI = 

2=1 

i.e., it is odd when exactly one ji is odd. We start from the constant function 

A (0) (a;) = l (5) 

for all x, where 0 = (0,..., 0). For definiteness we set X^\x) = 0 whenever ji < 0 for some i. 
Then we define the formal power X^ for admissible j with nonnegative indices in the following 
recursive manner: 

f IJI J i j odd, 

x^ = \ d _ ( 6 ) 

These are all functions on [x\,X 2 \- Note that when J is odd, the index i referred to in the 
first clause of (j6j) is unambiguously defined. The interdependencies of the X^ 1 are illustrated 
for d = 2 in Figure |T| We will say that the degree of X^ is |j|. 

To motivate to some extent what we have done we give the following relationship. 

Lemma 1 Let uo be a nonvanishing function on [ 11 , 2 : 2 ] and suppose that Luq = 0. Then for 
any nonnegative even multiindex 2 n, 

d 

L(uoX<W) = 2|n|(2|n| - 1) u 0 ^ nX^- 2 ^. 

i=l 

Proof. As a consequence of Luq = 0, is easily seen that the operator L admits the Polya 
factorization m 

L = — dpul d — 
uo u 0 

where d = d/dx and the functions in this expression refer to the corresponding multiplication 
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operators. Thus by ©, 


L{u 0 X (M) ) = —dpuld— ( 2|n|«o ^ 

uo Mo y J put i 

= 2|n|-^<9 j . 


-Si) 


A second application of © yields that this is equal to 

2|n|J- (^2\n\-l)J2nu 2 0 X^~ 2 ^ 

as required. □ 

The number cy of summands of the form © comprising XW is the same as the number 
of paths which advance (i.e. from predecessors to successors) from to and can be 

described as follows. We define cy = 0 when some j, < 0. Clearly eg = 1. Then by © we have 
recursively 

C-. 


yy_&) Jodd, 

d 

ES--SV J even. 


(7) 


i=1 


r2 » p > A' (0 ’ 2) r2 r- x(°’ 3 ) P > X(°> 4 ) r2 » X(°’ 5 ) 


x(b°) 


X(b2) 


'' 

X(b4) 


p 


X ( 2 ’°) r2 » x^ 2 ’ 1 ) v ' X (2 ’ 2) r 2 x( 2 ’ 3 ) P ^ X( 2 > 4 ) r2 > X( 2 ’ 5 ) 


X( 3 '°) 


X(3> 2 ) 


X( 3 > 4 ) 


X(4,0) r2 „ x(4.i) p ^ x<- 4 ’ 2) r2 r- X< 4 ’ 3 ) P ^ X( 4 > 4 ) r2 > x( 4 ’ 5 ) 


x( 5 '°) 


X( 5 > 2 ) 


x( 5 > 4 ) 


Figure 1: Construction of X^\ 
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( 8 ) 


By induction via predecessors it is readily seen that 


c-r = 



Consider a single nested integral appearing as a summand in The number of integra¬ 

tions following division by pu/j is [|J*|/2], while for each i (1 < * < d), the number of integrations 
which follow a multiplication by ?yug is easily seen to be [(ji + l)/2]. Here and always [a] means 
the least integer no greater than the real number a. One may verify that these indeed sum to 
|J|. Define 

M 0 = sup - — 21 , Mi = sup \riUo\. (9) 

[xi,x 2 ] \P U 0\ [xi,X 2 ] 


Lemma 2 The foi'mal powers satisfy the growth condition 

~,-n ri£il riiiil r j2+n [ Ja +1 1 

\X^\x)\<c t M ^ 2 l M[ 2 J M 2 L 2 • • • Mf J 


X — Xq\^ 


for X\ < x < X 2 - 


( 10 ) 


Proof. Suppose that j is even. Then by the inductive hypothesis 

\X^\x)\ < \j\ j {sup-\-)J2\X^- di) {t)\dt 

Jxo \P U 0\ i=1 

rL i T 3d + 1 1 \ 

\j\\t — u?o | 1—1 dt. 

We then integrate and note that \ji/2\ = [(ji + l)/2] since all ji are even, obtaining the bound 



r x ( * ri£Hii [ii+ii 

[lil 

[^1\ 

< 

/ 2 M 2 

• • m\ 2 J • 

j 


Jx o \i=1 


/ 


] 71df 2+ ^ ^ 


xo\ j+k 


\i=l 


which by reduces to (flUl) . The verification for \j\ odd is similar and indeed simpler. 


1.2 Construction of 

The construction of X^ 1 in terms of nested integrals of the form 0 is analogous to that of 
X^\ However, there are some notational complications. The indices could be handled in 
various ways; our choice, perhaps purist, is as follows. Now the j will have entries with common 
fractional part ji — [ji] = 1/d. We will call j admissible when at most one of the the integral 
parts [ ji\ is odd, while the parity of j is again that of the integer | j\. 

To start the recursion we use the predecessors of (l/d)l , i.e. j = (1/d, 1/d,..., 1/d, — 1 + 
1/d, 1/d,..., 1 /d), defining the constant functions 

X «i/d)i-gi), x s = I (H) 

a 
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for i = 1 These are functions of degree 0. For definiteness we define X^(x) = 0 

whenever some ji < 0 except as specified by m ■ The formal powers for the remaining 
admissible j > 0 are defined by 


X&> 



j even, 


J odd, 


( 12 ) 


as outlined in Figured] This formula differs from © not only in the exchange of even and odd, 
but also in that the indices and coefficients have different interpretations. One justification for 
this notation is the role of the degree |/|, cf. Lemma [I] 

Analogously to Lemma [1] we find 


Lemma 3 


L(u 0 X( 2H +^) = (2(|n| + l)(2|n|)n 0 ^r i X( 2 "- 2?i +^). 

i= 1 


In verifying this it is useful to note that |(l/ci)l| = 1 and to use the linearity of the degree 
operator | • |. 

The number cy of terms in X^) is determined recursively by setting = 1, while 

otherwise cy = 0 if some j t < 0, and then 


J even, 

d 

^ odd ’ 


(13) 


which is analogous © but again with a different interpretation of the indices. The number of 
integrations following division by puy is now [(|J*| + 1)/2], and those following multiplication by 
riUQ number [(ji — 1 /d + l)/2]. This gives the growth estimate: 


Lemma 4 


f j +1 1 f ji — 1/d+l 1 f 3 d.- 1 / d + 1 ] 

\xW(x)\<c f M [ 0 2 2 ■ ■ ■ 2 J 


X — Xq 


m 


6 






X ( 2’-2> ->. A ( 2’2> 

P 




1 i) r 2 


. xd’ 1 ^ P >- X<i’ 2 i> r2 A ( i’ 3 i> p A'(i’ 4 i } 


A (1 2 > 2 2 ) 


A (1 5’ 4 5) 


X< 2 5’3) r2 > x^i’ 1 ^ P > X (2 i’ 2 i) r2 > A' (2 i’ 3 ^ P > A< 4 i’ 4 i> 


A< 3 i-i> 


A (3 7’ 2 ^ 


A^ 3 5’ 4 3) 


A( 4 i’i> r2 > a( 4 5' 4 5) p > A( 4 i- 2 3> r2 > A^ 4 5’ 3 i) p > A( 4 i’ 4 i> 


X^i’i) 


X <5 b 2 i) 


A (5 5' 4 5> 


Figure 2: Formal powers X *^1 for d = 2. 

2 SPPS series and characteristic function 


2.1 General solution 


We define the SPPS functions U \, U 2 as 


Ml 


u 2 


u 0 


E ' y(2n) \ rai \n d 

^ir?ni A Al ’" Xd ’ 

x (2n+^l)yii 


n>0 

u 0 J2 

n>0 


(2|n|)!' 

1 

(2|n| + l)!' 


(14) 


where the sums are over all nonnegative multiindices n. Note that the degree of A'( 2n+ ^ 1 ) is 
2\n\ + 1. The main result, which mirrors that of [22], is as follows. 


Theorem 5 Let p,q,n,... ,ra be continuous complex-valued functions of the real variable x £ 
pro, Xi], with p continuously differentiable and p(x) 0. Let the differential operator L be 
defined by @1. Then for every A = (Ai,..., Ad) £ C d the two series in 0 converge uniformly 
on x £ [xq, xi], and the functions ui, U 2 thus defined are linearly independent solutions of {Ip. 
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Further, their derivatives are given by uniformly convergent power series, 


A = ^v + —v 1 

1 Uq 7mn ^ (2177.1 — 1 )! 1 Z 


pu0 (2 I”I ^ 1)1 i=l 


4 = ^w 2 + — V— 

«o pu 0 ^ 0 (2|n|)! 1 d 


( 15 ) 


For every va/ue o/ A, t/ie initial values of these functions are equal to 

ui(x 0 ) = u 0 (x 0 ), u'^xq) = zto(x 0 ). 


w 2 (a:o) = 0, u' 2 (i r 0 ) = 


p(a’ 0 )Mo(a;o)' 


( 16 ) 


Proof. This proof is quite analogous to the proof for d = 1 given in [23], but certain details must 
be taken into account when there are more spectral parameters. First we verify the convergence. 

Let A = max(|Ai|,..., |Ad|). Recalling (O, let M = max(M 0 , Mi ,..., Mu). Now by (fTUll . 
\X (2 ^{x)\ < c 2 \n\M [ W ] ■ A/[ ni ]+-+M \ X2 - Xl | 2 l"l = CfMW\x 2 - 2q| 2|s| 
so by ([H]) the summands in the formula for u\ in (TUI) are bounded by a 2 l™l/(2|n|)!, where 

a = \JMA\x 2 — x\\. 

Since a 2 \ n \ = a 2ni a 2n2 ■ ■ ■ a 2nd , we can factor the sum Y^’o’ a /(2|f?|)! into a product of d sums, 
each of which is equal to cosh a. By comparison with this finite sum it follows that the series 
for ui converges uniformly to a function bounded by coslr d a. By similar arguments the series 
for u\ , « 2 , and u' 2 also converge uniformly, and this justifies the term by term differentiation. 

By Lemma [T] 


L m = Y 


| n |=0 


(2|n|)! 


A” 1 ■ ■ ■ \ n d d L^uoX^) 


°° \n i \n d d 

E - 1 j5 H jF(2|«l)(2|S|-l)»oE(r.xP 


n— 28 i 


| n |=0 


°° \rai \ n d 

Un V 1 "' Ad V r .y< 2 "- 2<5 b 

U ° L. (2|n|-2)! ^ UX 

| n |=0 V 1 1 ' *=1 


Rearrange the last double sum as 

d 


a °° \ni \n d 

Vr V 1 — d y(2(a-ft)) 

^ 1 ^ (2(|n|-l))! 


2=1 | n |=0 
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and reindex each nt down by 1 , using the assumption that 


0 when j, : < 0 for some i: 


E 


l«l=o 


\rii \rid _ 

A 1 y(2(n—6i)) 

(2(|n|-l))! 


E 

| n |=0 


A ? 1 


■X 


m +1 


••A" 


( 2 (|n|))! 


X ( 2 n) 


Thus we have 



| n |—0 


\™i \rid 

A 1 "~ A d -y-{2n) 

(2(|n|))' 


Lui = 



u 1 


and the same argument verifies the corresponding statement for 112 ■ The final statement re¬ 
garding the initial values follows from the fact that X^\xo) = 0 whenever even a single ji is 
positive, so only the constant terms survive in the series for u \, u [; similarly all but the lowest 
degree terms involving X^(xq) also vanish. □ 


It is well known that when p , < 7 , Tj are real-valued, a complex nonvanishing solution uq of 
Ly = 0 can be obtained as a complex linear combination of any two linearly independent 
solutions; in fact, by considerations of dimension one sees it is not necessary for the coefficients 
be real-valued for such a nonvanishing solution to exist. The hypotheses of Theorem [5] could 
be weakened, for instance by only requiring 1 /{pvffi and to be continuous, but we will not 
enter into such details here. 


Corollary 6 Let 111,112 be given by G2P- Define 


Vl 

V2 


1 

u 0 (x 0 ) 


Ml - p(x 0 )u' 0 (x 0 ) u 2 , 


p(x 0 )u 0 (x 0 ) u 2 - 


Then v\,v 2 satisfy the normalizations 

vi(x 0 ) = l, v[(x o ) = 0 1 

v 2 {x 0 ) = 0, v' 2 (x 0 ) = l. 

Observe that v±,V 2 are also represented as power series in Ai,..., A<p 


2.2 Generalized Sturm-Liouville equation 

We consider now equations of the form 


d 

Ly = E A iRi[y\ 

1=1 


where we define 

Ri[y] = fiy + Siy 1 


(17) 


(18) 
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for given functions r ? ;, Sj, i = 1 Thus © is the particular case where all s,; vanish 

identically. In [26] SPPS formulas were developed for OB for the case d = 1. The multiparam¬ 
eter version is as follows. The formal powers are taken now with the same starting values as 
previously but with the followng modified recursive definition: 


x& = 1/1 

1 u 0 Ri[u 0 X^ E (J odd); 

X™ = I/I J 

f 

P<ti 

s 

II 

El 

^ UoRiluoX^-^} (j even); 

* 

s 

II 

El 

( 

p u lU 

again we use integral entries in J* for X^ and 

non-integral entries for X^h 


We will need the common bound 

1 


M = sup 


, M 0 f?i[uo] , •••, \u 0 Rd[uo]\, 


[x7j 2 ] \ \ p u o i 

of the functions appearing in the considerations below. For odd j it follows from (fT51) (fl9l) that 


RiluoXV-^} = R^X^ 


-Si) , (1/ I ~ 1 ) s » 


pu 0 


E* (f 


-Si-S,,) 


( 20 ) 


which implies that the formal powers may be calculated without recourse to numerical 
differentiation (other than for u o) and that 


/ d 

u 0 R t [u 0 X ( ^ } ] 

v i '=1 ' 

Analogous statements hold for X^K 

Lemma 7 For all x G [x 2 ,X 2 ], the inequalities |X^| < P|j| (x) and IA'^1 < P\f\{ x) hold, 
where 


< M \XV -**| + (|J| - 1) E I • (21) 



zo 


Xo\ 


for integral j > 0. 


Proof. First we consider X^\ i.e. j has integer entries. Write Ek = (M \x — Xq\) k /k\ so 
| M f Ek- 1 | = Ek . The inequalities are clearly valid when |/| is 0 or 1. Suppose that it is valid 
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for |/| up to n — 1. Now if |/| = n is odd and / has an odd entry in the *-th position, we calculate 
that 


n — 3 


P n -i{x) = d^(n-l)! V ( 2 W 

, n +1 \n-l- kj 

k =-T- 

n —2 ✓ 

d(n - l)|P n _ 2 (x)) = d“(n-l)! ^ f 

fc_n-l V 


n — 3 


,-2-k 


Then by the inductive hypothesis and (Uhl) . (12T1) . 


X^(x) 


< n 


[ M(P n -i(x) + d(n- l)P n - 2 (x)) dx 

J Xq 

/ n —2 

d 2 n\ | E n+i + 


n — 1 


n—3 


— 1 — k. 


n — 3 


,-2-k 


Ek+i + E n 


n —1 / _ _ \ 


= P n {x) 


as is seen after reindexing k + 1 to k and then noting that [n/ 2] + 1 = (n + l)/2. On the other 
hand, if / is even, then a similar, simpler argument verifies the inequality. 

The verification for X ) is analogous. □ 

The following results for the generalized formal powers are now proved in exactly the same 
way as Lemmas [1] and [3] and Theorem [5] 


Lemma 8 


and 


L[u 0 X (2n ' 


L(u 0 A'( M+ 3 r )) 


d 

= 2|n|(2|n|-l)^i? i [ Uo X( M - 2 ^)] 
2=1 

d 

(2(|n| + l)(2|n|) 

2 = 1 


Theorem 9 Let p, q, r ±,..., r^, si,..., Sd be continuous on [xo, Xi], with p continuously differ¬ 
entiable and p{x) 0. Define X® and X^P by 1191) . and then define u\, U 2 by These 

series converge uniformly on x £ [xo, Xi] for every fixed A = (Ai ..., Xf) £ C d , and are linearly 
independent solutions of the generalized Sturm-Liouville equation m- Their derivatives are 
given by fU and they satisfy the initial conditions ([ZIP- 
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2.3 Spectral problems 


The treatment of multiparameter spectral problems by the SPPS approach is the same as for a 
single spectral variable. Consider for simplicity linear boundary conditions of the form 

au(ari) + aV(a:i) = 0, Pv(x 2 ) + P'v'{x 2 ) = 0. ( 22 ) 

For the general solution v = C\V\ + c 2 V 2 with vi, V 2 given by Corollary [ 6 l this gives rise to a 
system of two equations in C \, C 2 with determinant 

a(/3v 1 (x 2 ) + /3v' 1 (x 2 ) - a'ffiv 2 (x 2 ) + /3v' 2 (x 2 ). 

Thus (1221) is satisfied when y(A) = 0, where 

x(A) = — a'/3 vi(x 2 ) + ct/3 v[(x 2 ) — a '/?' v[(x 2 ) + a/3' v' 2 {x 2 )- (23) 

Theorem [5] represents y(A) as a power series in A = (Ai,..., A d). Solutions of the boundary 
value problem are precisely the zeroes of this analytic function of several complex variables. 

In like manner, nonlinear or mixed boundary conditions will also produce a characteristic 
function. When these conditions are analytic, the result will be expressible as a power series in 
the A i, although it may be more convenient to leave it as a function defined as a combination 
of power series with other types of functions (cf. (l28l) below). 

Similarly, one may impose boundary conditions at more than two points. One way of solving 
such a problem is by converting it to an integral equation Huang. With the approach described 
here, one simply evaluates the SPPS representation at all boundary points required, in order to 
obtain the desired set of simultaneous characteristic equations. 

2.4 Remarks 

2.4.1 Reduction to simple cases 

We note that for d = 1 (i.e. (0) = (0), (1) = (1)), the starting integral of the X ^ family is 
X <0> = 1, and by (fill) for the family X^ 1 the starting integral also reduces to 

= X (0) = - = 1 

1 

Further, for general d the degree-1 power is simply the integral 



which coincides with in the case d = 1. Thus our notation is consistent with the “classical” 
definition of [23] . 

Considering d > 1, let us suppose that r,; is identically zero for every i / ip- Then X^ 
will vanish whenever j contains a ji > 0 where i ^ iq. The surviving powers 
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form the sequence X^' 1 of classical 1-spectral-parameter formal powers in the single variable 
A t. Similarly, the X^> reduce to the sequence and the series u\, 112 become the classical 

SPPS solutions. 

On the other hand, when all the r; are equal, the formal power X^ 1 is unchanged^when 
the entries j 1 ,... - jd are permuted, so the sum only depends on the degree |J*|, giving I® = 
and similarly = cjX^^\ which are multiples of the classical formal powers. It 
follows that Hi, U 2 are the classical solutions obtained using |A| = Ai + ■ ■ ■ + Ad in place of the 
single spectral parameter. 


2.4.2 Computational aspects 


We make a few observations to simplify the task of programming the formal powers. One can 
omit the factors |J*| in the recursive definitions (0, (fT^l) of and X^\ producing “rescaled 
* * 

powers” X ^ and X^> defined by 


X& = 



3 odd, 


j even. 


* 

and similarly for X^\ Then by induction 


X& = — -X (J) 

|J1- 


X& = —X (f) . 

|J 1 - 


Besides this saving in multiplications when calculating the formal powers (and often avoiding 
calculating with very large numbers), it is no longer necessary to divide by these factorials to 
obtain the terms in the sums for u \, 1 x 2 , u'i , u 2 ; i.e., we have simply 


u 1 =uoY, x( - H) XT ■■■*?, 

n 

etc. This is because the coefHcent of each formal power in the formulas (fl4l) is precisely the 
reciprocal of the factorial of its degree. 

The construction of the tables for X and X is seen to be identical when we disregard the 
initial terms X^ 1 ^ d ^ 1 ~ Si ' ) from the second table. That is, according to whether we insert the 
function 1 = or f l/(pug) = X^ 1 ^ d ' >1 ' > in the upper left hand corner, the same procedure 
of multiplying and then integrating will produce the entire table for X^ or X ^ respectively. 
Both tables and the corresponding power series can thus be computed via a single program, 
except that in the formula (1151) for u 2 , the first term corresponding to ft = 0 contains negative 
exponents and is not found in the truncated table. Its value is 


1 

pu 0 


0 ! 


^x ((1/d)r 


Si) 


\0 

A 1 A d 


2 = 1 


1 

pu 0 ' 
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This term must be added in separately to obtain u 2 . 

When programming, one may likely prefer to drop the fractional parts of the indices, using 
effectively 

XiT) = jfOr+CVdlT). 

In the development of the theory given above, this amounts to replacing the coefficient |J*| with 
|J| + 1 , which is the true degree of X^). 

It is easily seen that if uq is a solution of (JT]) for a fixed multiparameter (Ai,o, • ■ • , A^o) £ C d , 
then our construction of X^\ X ) w ill produce series in powers of Ai — Aqo, ■ • ■, A d ~ A^.o 
analogous to m (HU- This can be used to recenter the series for obtaining increased accuracy 
as in [25] . 

When calculating one must truncate the problem, say by using a finite number M of points of 
[x\. X 2 ] when integrating, and by approximating the series (1141) (1151) with polynomials formed of 
the terms for |n| < N. The total number of formal powers in {X^, X ^}^<jv grows as 0(N d ), 
so the memory requirement is of the order 0(MN d ). For boundary value problems this can be 
reduced by saving only the last value X^(x 2 ), X^\x 2 ) once the values interior to the interval 
are no longer needed for further integrations. The resulting memory cost O(MN) + 0(N d ) is 
in fact a great savings since often M is much larger than N. 


3 Numerical examples 

We give some examples for d = 2. The operational parameters M, N are as described at the 
end of the last section; calculations were carried out in Mathematica. 

3.1 Boundary value problems 

Example 1. This simple example uses constant coefficients p = 1, q = 0, rq = r 2 = — 1. 
The equation u" = — (Ai + A 2 )u has normalized solutions iq(:r) = cos(-\/Ai + \ 2 x), V 2 {x) = 
sin(\/Ai + \- 2 x)/\/\i + X 2 - On the interval [x\,X 2 ] = [0, 7 r], the SPPS solutions of Corollary [ 6 ] 
with M = 800, N = 20 are found to agree with these formulas to within 10 -9 for |Ai| < 1. 
As is common with polynomial approximations, the accuracy drops rapidly for larger values 
of |Ai| when the truncation limit N is fixed. We impose the boundary conditions u(0) = 0, 
u( 7 r) = 0. The graph of the characteristic function y(Ai,A 2 ) (eigensurface) is shown in Figure 
[3 The eigencurves % = 0, calculated numerically from % via the function ContourPlot in the 
figure, coincide with the solutions of 


Ai + A 2 = 


k 2 n 2 

b 2 


for k = 1,2,3,4. Indeed, the values of |x(Ai, A 2 )| for |Ai| < 5 for k = 1,2 ,3,4 are less than 
10 -12 , 10~ 12 , 10 _1 °, 10 -5 respectively. When the maximal degree of the powers is reduced to 
N = 16, the level curve for k = 4 is visibly far off the mark. 
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Figure 3: Characteristic function and zero-level curves (Example 1). 


Example 2. This example, with p{x) = 1, q{x) = cos a;, r\(x) = cos(ar), r 2 ,{x) = cosx, 
which is not amenable to a solution in closed form, is chosen to illustrate level sets which are 
not connected and which contain closed curves. Using the same interval [0,7r] and boundary 
conditions u{0 ) = 0, u{n) = 0, we find the characteristic function and its zero set as depicted 
in Figure [4] For illustration we take an arbitrary section A 2 = 1.0, and restrict x t° this 
value (Figure[5|). The corresponding numerical pairs (Ai, A 2 ) determine an ordinary differential 
equation which can be solved numerically by NDSolve using the boundary condition at x = 0 to 
define an initial condition. The resulting values at x = tt were found to differ from x(Ai, A 2 ) by 
less than 10 -6 when the experiment was carried out with M = 100, N = 12. The calculation 
of the characteristic function took about 0.3 seconds, and then each value of Ai less than a 
thousandth of a second on a portable computer (this does not include the time for checking 
by solving the initial value problem). The three eigenvalues Ai « —9.5644, —4.3944, 3.9177 in 
the range considered are easily located by techniques of numerical approximation of zeroes of 
polynomials. 

Example 3. The following example involves consideration of complex eigenvalues. The boundary 
value problem 

y"(t) + (E + zsgat)y(t) = 0, y(-l) = y(l) = 0, (24) 

where sgn.T is the sign of x, was studied in detail in m- a spectral surface is formed of pairs 




-10 -5 0 


Figure 4: Characteristic function and zero-level curves (Example 2). 
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(E,z) G C 2 . Let Ai = E, A 2 = z, ri(x) = —1, r 2 (x) = sgn(x) for — 1 < x < 1. 

For the SPPS calculation, due to the jump singularity in r 2 , it would be appropriate to 
integrate separately on [—1,0] and [0,1]. For this example, however, we simply calculate the 
formal powers with M = 10,000 mesh points, and settle for about five significant figures in 
the integrations. Since |y| does not change sign near its zeros, we take the logarithm; then the 
plotting routine (Plot3D) easily reveals the set where y(Ai, A 2 ) = 0 as shown in Figure[Gl where 
we have taken Ai real and A 2 purely imaginary. 

In [30] certain curves in the spectral Riemann surface were explicitly parametrized as 

Ai(s) = s 2 — h(s) 2 , A 2 (s) = 2 ish(s) (25) 

where s G U^LoK 71 S' 1/2)tt, (n + l) 7 r] and where h is defined implicitly by the relation 

s sin( 2 s) + h(s) sinh( 2 /i(s)) = 0 . 

These curves were used to show that the surface is connected by joining various 1-complex- 
dimensional parts. The first interval s G [tt/2, 7 r] corresponds approximately to 2.467 < Ai < 7 r 2 , 
0 < A 2 /i < 4.475, and is the smallest eigencurve revealed in the plot. For the values given by 
(1251) with values of s in this interval we find numerically that |y(Ai(s), A 2 (s))| < 10 4 for 
s G [7 t/2,7 r] when N = 20. 




Figure 6 : log |x(Ai, A 2 )| for (l24l) with N = 40 (left); detail of region around smallest eigencurve 
with N = 16 (right). 
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3.2 Application to electromagnetic transmission 


Example 4- This example is based on m from which we restate the minimum possible of 
background material. The plane R 2 = {(£,?/)} is partitioned into the regions 

Sli = {x < 0}, = {0 < x < &}, fi 2 = {x > b}, 


which are assumed to be composed of materials such that the index of refraction in f2i and fi 2 
takes constant values denoted ni, ro 2 respectively, while in the inhomogeneous region fig it is a 
function n = n(x) independently of y. These values are bounded below by 1. An electromagnetic 
wave of the form e _lfclX travelling in strikes the boundary line x = 0 with flo at an angle 9 
from the perpendicular, and is partially reflected back into Sli as u(x) = e~ lklX + Re lklX and 
partially transmitted into H 2 at x = d as u(x) = e~ xk2X . The parameter 

/3 = ksin9 (26) 

is introduced, where k = 2 tt/\ is the wave number in terms of the wavelength A (here A will 
not denote an eigenvalue). In flo the wave is governed by the differential equation 

u"(x) + ( k 2 n(x ) 2 — /3 2 )u(x) = 0 (27) 


(for the “s-polarization”, and a similar equation for the “p-polarization”). The problem is the 
determination of the complex constants R and T, known as the reflection and transmission 
coefficients. In m the formulas 

_ — k\k 2 v 2 {b) — u^(6) — ik 2 Vi(b) + ik\v' 2 (b) 

(v'i(b) — fcifc 2 u 2 (&)) + i(k 2 Vi(b) + k\v' 2 {b)) ’ 

= 2ik 1 {v 1 {b)v' 2 (b) - v , 1 (b)v 2 (b))e- ik2b 

{v[{b) - k\k 2 V 2 {b)) + i{k 2 vi(b) + kiv' 2 (b)) ’ 


were given, where Aq = y/k 2 nf — /3 2 , k 2 = y/k 2 n 2 — /3 2 . It was shown how by fixing k in (l26l) 
and then using /3 2 as the spectral parameter, the SPPS formulas for dimension d = 1 can be used 
to calculate R and T for varying angles of incidence 9. Examples were given for three sample 
functions n(x). All were for normal incidence j3 = 0, for which it is not difficult to calculate the 
solution of the differential equation analytically in terms of special functions for the examples 
considered (see for example [33]), and thus compare the accuracy. Similar calculations using 
SPPS were carried out in |12| . again for normal incidence, with many graphs comparing the 
results to other numerical methods used in optics. 


Equation ([T]) for d = 2 with Ai = /? 2 , A 2 = —k 2 , r±(x) = 1, r 2 (x) = n 2 takes the form (l27l) . 
We apply Corollary [5] to obtain normalized solutions vi,v 2 , and then substitute these together 
with _ _ 


ki 


—Ai — A 2 n 2 , k 2 = J —Ai — X 2 n 2 


(29) 


in ([28]) . This produces analytic expressions i?(Ai, A 2 ), T{ Ai, A 2 ) which, while they are not simple 
power series, serve conveniently for calculations. 


We will take one example, the “hyperbolic” refractive profile 


n{x)=n{ O)e (a:/b)los(n(6)/n(0)) 


(30) 
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Figure 7: Logarithmic plots of |i ?| 2 (solid), (ri 2 /ni)|T | 2 (dotted) and formula (IdTl) (dashed) as 
functions of the adimensional magnitude 6 /A 2 € [(3/(2tt), 1 ], for n{x) given by (1301) . 


with d = 1, 7 ii = 1.0, n(0) = 1.4, n[b) = 2.1, n -2 = 1.5. Further, we set b = 1. In Figure |7] all 
graphs were plotted after a single calculation of the series for x(Ai, A 2 ) and its substitution in the 
expressions (l28l) for the parameters given above. It follows from (l26l) that 6 /A > /36/(27t), which 
determines our starting point for plotting the curves. For normal incidence (3 = 0, conservation 
laws require the expression 

|7 ?| 2 + —|T | 2 (31) 

ni 

to be equal to 1; this is seen in the first graph, which agrees with Figure 6 of [TO]. For other 
values of [3 we have spot-checked numerically by selecting various values of the dimensionless 
quantities (3 and 6 /A, then solving the corresponding (1271) numerically with NDSolve, as in the 
previous example. The final values v\(b),... , v' 2 {b) produce values of R, T via (l28l) for checking 
against the x-values plotted here. The results are given in Table [I] All of the data here is 
affected by the fact, observed in 10 }, that arithmetic operations in (128)) reduce the accuracy 
produced by the differential equations by several significant figures. 


4 Closing remarks 

We have shown how the representation of the solutions of the Sturm-Liouville differential equa¬ 
tion in terms of power series in a single spectral parameter may be generalized to several 
parameters Ai,..., A^. We hope that this will make possible a deeper analysis and simplified 
computation for many problems in physics and engineering, which have been approached up to 
now by fixing the values of all parameters but one, and solving by uniparameter methods. 

Regarding the many aspects of uniparameter SPPS theory which have been developed up 


18 

































M = 30, N = 10 M = 50, N = 16 


p 

b/X 

0.01 

0.05 

0.1 

0.15 

0.2 

P 

b/X 

0.01 

0.05 

0.1 

0.15 

0.2 

0.005 
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6 

7 

4 

0 

0.005 

7 

7 

7 

7 

3 

0.01 

— 

6 

7 

4 

0 

0.01 

— 

7 

7 

7 

3 

0.1 

— 

— 

7 

4 

0 

0.1 

— 

— 

7 

7 

3 

0.5 

— 

— 

7 

4 

0 

0.5 

— 

— 

7 

7 

3 

1 

— 

— 

7 

4 

0 

1 

— 

— 

7 

7 

3 


Table 1: Number of significant digits in y(Ai, A 2 ) for selected values of /3, b/X. 


to now, we point out as illustrative examples only two possible areas for using several spectral 
parameters. 

The so-called Sturm-Liouville pencils 

{pu'Y + qu= £>A*)« 

have been investigated from the SPPS perspective in arXiv: 1401.1520. This equation is a 
particular case of JU with Ai = 1, A 2 = A, ..., Ad = A d_1 . Thus our formulas provide the SPPS 
series for this equation directly. 

In another direction, coefficient functions with singularities at one of the endpoints [x\, X 2 ], 
such as occur in Bessel’s equation, have led to modified versions of the SPPS formulas [15- 
Similar results can be expected to hold also for several spectral parameters. 

We close with the observation that an alternative construction to the one described in this 
paper may be developed by first setting all but one of the spectral parameters to zero, for 
example considering 

( py'Y + qy = A my, 

and writing down the classical formulas for solutions ui[ Al ^, depending on this parameter. 
These can be regarded as solutions of 

( py'Y + (9-Aid )y = 0, 

and after choosing a suitable nonvanishing linear combination, this can be used as the seed for- 
solving 

{py'Y + (<? - A in)y = A 2 r 2 u 

to obtain w[ Al,A2 ^, and so forth. Even for the case d = 2 the resulting calculations to 

recover the coefficients of the SPPS series turn out to be surprisingly complicated, and involve 
many products of the nested integrals which cancel out at the end. The author is grateful to S. 
Torba for suggesting the simpler approach followed in the present work. 

This research was partially supported by grant 166183 of CONACyT (Consejo Nacional de 
Ciencia y Tecnologfa), Mexico. 
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